Numerical modeling of 1-D transient 
poroelastic waves in the low-frequency range 



Guillaume Chiavassa , Bruno Lombard 6 , Joel Piraux & 

a Ecole Centrale de Marseille and MSNM-GP, 13451 Marseille, France 
b Laboratoire de Mecanique et d'Acoustique, 13402 Marseille, France 



Abstract 

Propagation of transient mechanical waves in porous media is numerically investi- 
gated in ID. The framework is the linear Biot's model with frequency-independant 
coefficients. The coexistence of a propagating fast wave and a diffusive slow wave 
makes numerical modeling tricky. A method combining three numerical tools is 
proposed: a fourth-order ADER scheme with time-splitting to deal with the time- 
marching, a space-time mesh refinement to account for the small-scale evolution of 
the slow wave, and an interface method to enforce the jump conditions at interfaces. 
Comparisons with analytical solutions confirm the validity of this approach. 
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1 Introduction 



The propagation of mechanical waves in porous media is of interest in many 
areas in applied mechanics, including industrial foams, spongious bones and 
petroleum rocks. The most-widely-used model describing the evolution of 
small mechanical perturbations in a saturated porous medium is that proposed 
by Biot in 1956. Two regimes were distinguished by Biot: one corresponding 
to a low-frequency range [3], and one to a high-frequency range, where some 
of the physical parameters depend on the frequency pE]. We focus on transient 
mechanical waves whose frequency content lies in the low-frequency range. 
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Until the mid 90's, Biot's equations were mainly studied in the harmonic 
regime. Various time-domain methods have been proposed since, based on 
finite-differences pi23] . finite-elements jl0 ; 26j, boundary-elements [22], and 
spectral methods [B] . Since non- realistic values of the physical parameters were 
used, the real difficulties arising when performing time-domain simulations 
were often overlooked [7][TT] . These difficulties are induced by the coexistence 
of two solutions with radically different dynamics: a propagating "fast wave" 
and a diffusive "slow wave" |8J. The latter is highly dispersed and attenuated, 
and remains localized near sources and interfaces. 

The aim of the present study is to develop an efficient numerical method to 
compute the solution made up of these two waves. A time-splitting is used 
together with a fourth-order ADER scheme [TH] . A flux- conserving space-time 
mesh refinement is implemented at the places where the diffusive wave is 
localized [2]. Lastly, an immersed interface method gives a subcell resolution 
of the interfaces and accurately enforces the jump conditions between various 
materials [21]. The numerical tools used are described and tested in ID. 

The paper is organized as follows. The Biot's model is briefly recalled in section 
2. The numerical tools are described in section 3. Section 4 presents numerical 
experiments confirming the validity and efficiency of this approach. In section 
5, conclusions are drawn and some future perspectives are suggested. 



2 Problem statement 

2.1 Biot's model in the low-frequency range 

Biot's model describes the propagation of mechanical waves in a porous medium 
consisting of a solid matrix saturated with fluids circulating freely through the 
pores. The underlying hypotheses in the low-frequency range are as follows: 

• the wavelength of the perturbations is large in comparison with the diameter 
of the pores, as well as with the representative macroscopic volumes; 

• the amplitudes of the perturbations are small; 

• the elastic and isotropic matrix is fully saturated by a single fluid phase. 

The model is based on 10 physical parameters: the density pf and the dynamic 
viscosity r\ of the fluid; the density p s and the shear modulus p of the solid; 
the porosity < < 1, the tortuosity a > 1, the absolute permeability k, the 
Lame coefficient A/, and the two Biot coefficients (3 and m of the saturated 
matrix. The following notation is introduced: p w = ^ pf, p = (f) pj + (1 — <p) p s , 
and x — P Pw — pj > 0. The unknowns are the elastic velocity v s , the elastic 
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stress a, the filtration velocity w = 4>{ v f ~ v s), where Vf is the fluid velocity, 
and the acoustic pressure p. When dealing with heterogeneous media, all of the 
parameters are assumed to be piecewise constant and discontinuous across the 
interfaces. Denoting x = a the location of the interface and k s its hydraulic 
permeability, the jump conditions are [12]: 

[v s (a, t)] = 0, [w(a, t)] = 0, [cr(a, t)] = 0, [p(a, t)] = w(a~, t). 

(1) 

Since w is continuous, no privileged orientation is induced by the last condition 
in (pP). If « s — > +oo, the open-pore conditions are obtained, with perfect 
hydraulic contact [5]. If k s — > 0, no motion of the fluid relative to the matrix 
occurs: the pores are closed. 

In the case of a harmonic wave of frequency /, the low-frequency Biot's theory 
is valid as long as / < 0.15 f c |[4], with 



C 2 7T a K Pf 

When / > 0.15 f c , the width of the viscous boundary layer is smaller than that 
of pores, and the flow of fluid is no more of Poiseuille's type. Consequently, a 
non-polynomial frequential correction is required in the dynamic permeability 
k/t], leading to fractional derivatives in the time-domain [20J. 



2.2 Initial boundary-value problem 



Setting U = (v s , w, a, p) T , 

1 o o - Pw /x-Pf/x\ 

Pf/X p/x 

(A/ + 2/i) -(3m 
(3 m m 



V 



Tj_p 



^0-p//p0 0^ 

10 







(3) 

the Biot's evolution equations [3] are written together with the jump condi- 
tions ([1]) in the form of a first-order linear system with a source term 



at ox 
[CU(a, t)]=0, 

U(x, 0) = U (x), 



-SU ifx^a,t>0, 



(4) 
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where the 4x4 matrices C are derived from ([I]). A detailed mathematical 
analysis of the solution to (j4l) in a homogeneous medium can be found in 
[TO] . Here we restrict the discussion to basic properties that are directly useful 
for numerical modeling purposes. The spectral radius of S is R(S) = ^ ^. 
The eigenvalues of A are real and denoted ±ci and ±c 2 , with c\ > c 2 > 0. 
Injecting a mode ^ l - kx ) in (j^j), where u is the angular frequency and k the 
wavenumber, gives the dispersion relation 

Ak 4 + B(u) k 2 + C{u) = 0, 

A — k m (A f + 2 u — /5 2 ml , 

(5) 

S(cu) = — k ((A/ + 2/i) p w + m (p - 2(3 p f )) J 2 + i 77 (A/ + 2 /i) u, 

C{u) = x K ^ — i V P^ 3 , 

with roots ±ki and ±fc 2 , 3fJe{fci >2 } > 0. The phase velocities c\ = uj / 3?e{A;i} 
and c 2 = o;/3?e{A; 2 } are strictly increasing functions of frequency, tending 
in the high-frequency limit towards c\ and c 2 . The corresponding waves are 
called the fast wave and the slow wave, respectively. If 77 = 0, which is phys- 
ically irrelevant with usual media, the waves are purely propagated and the 
mechanical energy is constant. If 77 7^ 0, the fast wave is almost non-dispersive 
and non- diffusive. On the contrary, the slow wave becomes highly dispersive 
and diffusive. If / <C f c , then c 2 c 2 : the slow wave tends towards a non- 
propagating mode [8] . The direct contribution of the slow wave to the overall 
wave propagation processes is therefore negligible. However, the accurate com- 
putation of the fast wave depends crucially on the effects of the slow wave on 
the balance equations at interfaces [5]. 



3 Numerical tools 



3.1 Numerical scheme 



This subsection deals with the numerical resolution of (J4]) far from a, that is 
where the stencil does not cross the interface. A uniform grid is considered 
here, with the spatial mesh size Ax and the time step At. An approximation 
E/" of U(xi = i A x, t n = to At) is sought. The numerical methods recalled in 
section [1] usually consist of solving simultaneously the propagating part and 
the source term in (j3J). If rj 7^ 0, a Von- Neumann analysis of stability typically 
yields 

Ai£min (^'^))' (6) 
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which is highly restrictive since R{S) may be large. For instance, CFL = 
ci At/ Ax ~ 10~ 2 with the media considered in section (14.21) . and CFL « 
10~ 12 can be reached with highly dissipative fluids such as bitumen. A much 
more efficient approach is to split (J1J and to alternatively solve by Strang's 
splitting [TS] the propagating part 

§i u+ w- x AU = - < 7) 

and the source term part 

The equation ([7j) can be solved by applying any explicit two time step spatially- 
centered flux-conserving scheme. The latter can be written abstractly 

Ur 1/2 = H(UI S1 ..., U« +s ), (9) 

where s is the width of the stencil. In the present numerical experiments, a 
temporally and spatially fourth-order accurate ADER scheme was used [T9~|f2"3"] . 
with s — 2. This scheme accounts for waves over long distances with small 
dispersion and diffusion errors, even on coarse grids. The equation (jSJ) is solved 
exactly: p and a are unchanged, whereas the velocities become 



n+ l n+l/2 . Pf (-> n+1/2 n+ \ -11 £. £lL n +l/2 , ln x 

H — -[I — e K ^ 2 jw i , iu" = e K * 2 w . / _ ^Q) 



The splitting ©-(JH]) together with the exact integration (fit)]) yields the optimal 
stability condition CFL < 1. However, since A and S do not commute, it 
decreases the theoretical order of convergence from 4 to 2. In practice, the 
convergence rate measured was close to 3 in our experiments. 



3.2 Mesh refinement 



If rj 7^ and / f c , the slow wave has much smaller spatial scales of evolution 
than the wavelength of the fast wave. A very fine grid is consequently required 
to account for its evolution. Since the use of a fine uniform grid on the whole 
computational domain is out of reach in view of 2-D simulations, the grid 
refinement procedure provides a good alternative. In addition, since the slow 
wave remains localized near the sources and the interfaces, grid refinement is 
necessary only around these places. 

To limit the numerical dispersion on the coarse grid, it is preferable to also 
perform temporal refinement using a local CFL stability condition. Many al- 
gorithms for space-time mesh refinement have been recently developed in a 
variational framework, which ensure the stability of the coupling between grids 
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by conserving a discrete energy [16] . Here we adopt another approach based on 
flux conservation [2], which is more naturally coupled to the flux-conserving 
scheme (jUJ). The extrapolated values required to couple coarse and fine grids 
are obtained by performing linear interpolation in space and time on the nu- 
merical values at the surrounding nodes. In the case of the Lax-Wendroff 
scheme applied to the advection equation, the stability of the coupling was 
proved in [1] for arbitrary refinement factor. 

The mesh refinement with a factor q on a zone of width h introduces a com- 
putational extra-cost proportional to the product of grid nodes hq by the 
number of time substeps q. Consequently, one must estimate carefully q and h 
in terms of the physical parameters, in order to get efficient simulations. For 
that purpose, one uses the wavelengths Xx — ci/f and A2 = 02/ f deduced from 
(J5]). The numbers of grid nodes per wavelength of the slow wave, in the fine 
grid, and of the fast wave, on the coarse grid, must be roughly equal, which 
provides q ~ C\jci. This factor may be large: in section HJ one gets q ~ 64. 
Numerical experiments have shown that it is preferable to use intermediate 
smaller refinement factors, for instance 1 to 8 and then 8 to 64. Concerning 
h, the refined zone must contain a few wavelengths of the slow wave, leading 
typically to h = 4 A 2 in section |U 



3.3 Interface method 
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Fig. 1. Interface between media and Q\ and the irregular nodes around a. 

This subsection deals with the time-marching at the grid nodes where the 
stencil of flU]) crosses a. We define J by xj < a < x and we assume that 
the grid is uniform from £j_ s+ i to xj +s (figure [p. At the so-called irregular 
nodes xj- s+ \, ...,xj +s , the scheme ([9]) must not be applied naively, for three 
reasons. First, the spatial derivatives of the solution are not smooth across a, 
hence the order of convergence of the scheme decreases. Secondly, the subcell 
position of a inside the mesh is ignored, which leads to a 0( A x) error. Thirdly, 
the jump conditions ([1]) are not enforced in the scheme. The numerical waves 
diffracted by the interface may therefore not tend towards the exact waves. 

At the irregular nodes, we adapt an immersed interface method [21]. This 
method requires knowing the jump conditions satisfied by the spatial deriva- 
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tives of U . We deduce these conditions from 



d" 



dt r 



[CU(a, t)} 



d 



for all m > 1 
0, 



0" 



U)(a~,t), 



(11) 

where D m is a 4 (m + 1) x 4 (m + 1) matrix depending on the permeability k s 
and on the physical parameters around a. Let k be a positive integer. At the 
2 & grid nodes surrounding a, 2 fc-th order Taylor expansions of U(xi, t n ) on 
a^, together with the jump conditions (fTTI) . are written in the matrix form 



U(x J+k , t r 



( 



M 



U{a ,t n ) 



\ 



Q2k-1 
\ g x 2k- 



T U{a , t n ) j 



\0(Ax 2k )) 



'121 



where M is a 8 k x 8 k matrix. The Taylor series remainder term is removed 
from (I12p . and the exact values are replaced by numerical ones. The spatial 
derivatives estimated by performing inversion of (fl2l) are used to build smooth 
extensions of the solution, called modified values, on the right of a: 



J+1.....J + S, 



U* 



\2k-X 



.4, 



(2/c-i; 



M- 1 



( jjn \ 

u J-k+1 



V un j+k J 



(13) 

where J4 is the 4x4 identity matrix. Modified values on the left of a are defined 
similarly. They are then injected into the scheme at the irregular nodes: 



J — s + 1, J, 
J + l,..., J + s, 



UT 112 = Hn (Ut s , U n j, U* J+1 , W l+S ) , 



u 



n+1/2 



(14) 



Ho. U 



TT* TT n JJ n 

u Ji u J+l' •• ' u i+ 



where Hq q and Hq 1 denote the operator H with physical parameters of 
media f^o and respectively. Some comments about the interface method: 

• since the jump conditions are linear, the work is mainly carried out during 
a preprocessing step. At each time step, only small matrix-vector products 
(JTSJ) are required to compute the modified values. The computational cost 
is therefore negligible in comparison with that of time-marching; 

• the matrix M in (fl2|) depends on the jump conditions and on the position of 
a inside the mesh. Using the modified values (fl4l) introduces into the scheme 
a subcell resolution of the interface, which removes the O(Ax) error 
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• in the limit case where the parameters are the same on both sides of a and 
the hydraulic contact is perfect, U* = U™ if k > s. Consequently, we again 
obtain the scheme applied in homogeneous medium; 

• with a r-th order accurate scheme, the local truncation error of (I14p at 
irregular nodes is r-th order if 2 k — 1 > r [21]. The fourth-order ADER 
scheme therefore requires k = 3. As deduced from [13], k = 2 suffices to 
ensure fourth-order overall accuracy; 

• GKS analysis [H] has been performed in the case of inviscid fluids on a 
uniform grid to determine the stability of the hybrid scheme and (fT^|) . 
This analysis is based on the possible existence of discrete increasing modes 
emitted solely by the interface without any incident field [23] . A parametric 
study showed that the hybrid scheme is stable in tests 1 and 2 presented 
below, whatever the position of the interface, with k = 1, 2, 3. 



4 Numerical experiments 



Parameters 


^0 




Q 


n- L 


Pf (kg/m 3 ) 


1040 


1040 


1040 


10 


7] (Pa.s) 








10~ 3 


2.2 10~ 5 


Ps (kg/m 3 ) 


2650 


2211 


2650 


2650 


P (Pa) 


1.85 10 9 


3.54 10 9 


1.85 10 9 


1.85 10 9 


<t> 


0.3 


0.01 


0.3 


0.3 


a 


2 


2 


2 


2 


k (m 2 ) 


io- 12 


10 -16 


io- 12 


io- 12 


A/ (Pa) 


8.40 10 9 


4.69 10 9 


8.40 10 9 


2.43 10 9 




0.88 


0.01 


0.88 


0.35 


m (Pa) 


7.05 10 9 


2.46 10 11 


7.05 10 9 


5.37 10 7 


ci (m/s) 


2364.9 


2314.1 


2364.9 


1817.9 


c 2 (m/s) 


774.9 


1087.7 


774.9 


897.6 


ci (m/s) 


2364.9 


2314.1 


2364.5 


1817.8 


c 2 (m/s) 


774.9 


1087.7 


38.1 


30.3 


ci/c 2 


3.05 


2.12 


62.06 


60.0 


fc (Hz) 








22955 


52521 



Table 1 

Parameters and related data in tests 1-2 (central column) and 3-6 (right column). 
A 400-m domain is studied. Incident, reflected and transmitted waves are 
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denoted by I, R and T. Fast and slow waves are denoted by F and S. Except in 
test 6, analytical and numerical solutions are shown in solid lines and circles. 
Vertical dotted lines denote the position of mesh refinement. If 77 = in 
each media, computing the analytical solution is straightforward. Otherwise, 
it follows from a Fourier analysis. The source is a C 4 truncated sinusoid with 
a central frequency / = 30 Hz, as shown in Figure [2]-a. On the main grid, 
Ax = 1 m, and the computations are performed with CFL = 0.9. 

Two sets of physical parameters are used, shown in table HJ In tests 1 and 
2, they model sandstone saturated with water (Q ) and schist saturated with 
water (fl±), except that r] = [Hj. This is not physically realistic, but it sheds 
light on the ADER scheme and on the interface method: S = means that 
no splitting ([8]) is required. Nor is mesh refinement required, since the slow 
wave propagates. In tests 3 to 6, the parameters model sandstone saturated 
with water (Q ) and with gas Since t] 7^ and / <C f c , the slow wave is 
a static mode, which puts the focus on the mesh refinement. The refinement 
factor deduced from Ci/c 2 is 64: see section I3T21 and tabled) 



4-1 Inviscid media 



B 




100 150 200 250 

x(m) 




Fig. 2. Test 1: at initial instant (a) and after crossing the interface (b). 

First we consider identical media linked by an imperfect hydraulic contact 
k s = 10~ 16 m.s _1 .Pa _1 at a = 200.67 m. Wave conversions are shown in Figure 
O-b, and a good agreement is seen between the analytical and numerical values; 
on this scale, the reflected fast wave is not visible. Since no contrast of Biot's 
coefficients occurs, the unmodified scheme (Q without the interface method 
( TT4l) would propagate the incident wave without diffraction. 



Test 2 deals with media (Qq,, Qi) linked by a perfect hydraulic contact 1/ k s = 
0. Without using the interface method, GKS analysis and direct simulations 
show that the scheme is unstable and the solution grows exponentially when 
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100 150 



200 250 
x(m) 



(b) 




-3.5 -3 -2.5 

log(l/N) 



Fig. 3. Test 2: after crossing the interface (a) and convergence measurements (b). 

crossing x = a. Figure 0-a gives the analytical and numerical values obtained. 
The results of convergence studies are presented in Figure 0-b. As mentioned 
in section [3T3l k = 2 or k = 3 maintain the fourth-order accuracy of the ADER 
scheme; k = 1 does not suffice, leading to a rate of convergence of only 2.9. 




(b) 



: — » — » — & — 8 — e-i^^^m^^^ 






J 


K ' 










o o Q with refinement ? 






xxx without refinement 







192 196 200 204 

x(m) 



Fig. 4. Test 3: no refinement (a); zoom around x s , with and without refinement (b). 



4-2 Dissipative media 



Test 3 focuses on the homogeneous dissipative medium f2 excited by a ponc- 
tual stress source of finite duration at x s = 200 m. A direct discretization with 
no splitting would give CFL = 0.03 (jH]). A snapshot of p is shown in figure 
@]-a. Fast waves are advected rightwards and leftwards while the slow waves re- 
main localized around x s , and vary considerably on small spatial scales. Their 
numerical values are highly smeared. In figure |4]-b, the coarse grid solution is 
compared with a solution refined 64 times around x s . Good agreement is seen 
between the exact and refined numerical values. 
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Fig. 5. Test 4: snapshot of p at time t\ (a-c) and ti (b-d). Zoom around a (c-d). 

Test 4 is performed on (f) , fii), with 1/k s = 0. Two successive mesh refine- 
ments with factors 8 and 64 are used around a. Snapshots of p at t\ and t 2 > ti 
are shown in figure [5j During the interaction of the incident fast wave with 
the interface (a-c), since the slow waves have a greater amplitude than that of 
fast waves, they play a crucial role in the balance of momentum and mass: a 
poor assessment of the slow waves would invalidate that of the other ones. At 
ii (b-d), the slow waves are greatly attenuated and remain localized near a. 

Lastly, two examples of wave propagation across multiple interfaces are in- 
vestigated. The media Q and Q\ are repeated alternatively. Test 5 treats the 
case of two interfaces, with a known analytical solution (figure [6]). As in test 
4, one observes the small-scale evolution of the static slow waves generated by 
the interfaces. 

Test 6 deals with of 10 interfaces with randomly distributed position. 

Performing such a simulation has physical applications, even in ID [PT| . Figure 
[7] shows a at the initial instant (a) and after the wave has crossed the whole 
set of interfaces (b). For the sake of clarity, only grid refinements positions 
from 1 to 8 are shown. 
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140 1 60 180 200 220 240 260 280 195 200 205 210 215 

x (m) x (m) 



Fig. 6. Test 5: snapshot of p after crossing the two interfaces (a), zoom (b). 
(a) (b) 




50 100 150 200 250 300 350 400 50 1 00 1 50 200 250 300 350 

x (m) x (m) 



Fig. 7. Test 6: a at the initial instant (a) and after crossing the interfaces (b). 
5 Conclusion 

Numerical modeling of 1-D transient Biot's model was addressed here for waves 
whose frequency content lie in the low-frequency range. Three numerical tools 
were combined to obtain a method describing accurately the wave propagation: 
a fourth-order scheme with time-splitting; a space-time mesh refinement; and 
an immersed interface method. This method is required to account for the 
properties of the slow wave. Further research is suggested: 

• studying the high-frequency range [I], where k/t] is proportional to f 1 ^ 2 . 
Fractional derivatives are therefore involved in the time-domain [151120] : 

• accounting for dissipative effects in the solid skeleton; 

• coupling with nonlinear boundary conditions, to model seismic rupture; 

• extending the method to two-dimensional configurations. The validity of 
each of the numerical tools has already been established in 2D. 
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